{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b236b8d3",
   "metadata": {},
   "source": [
    "# Import an Extracted Network and Predict Transport Properties\n",
    "\n",
    "This example illustrates the process of both extracting a pore network from an image (that has already been binarized), then opening this image with OpenPNM to perform some simulations. The results of the simulations are compared to the known transport properties of the image and good agreement is obtained.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e13b67cf",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/amin/Code/OpenPNM/.venv/lib/python3.13/site-packages/porespy/filters/_lt_methods.py:21: FutureWarning: `square` is deprecated since version 0.25 and will be removed in version 0.27. Use `skimage.morphology.footprint_rectangle` instead.\n",
      "  strel = {2: {'min': disk(1), 'max': square(3)}, 3: {'min': ball(1), 'max': cube(3)}}\n",
      "/Users/amin/Code/OpenPNM/.venv/lib/python3.13/site-packages/porespy/filters/_lt_methods.py:21: FutureWarning: `cube` is deprecated since version 0.25 and will be removed in version 0.27. Use `skimage.morphology.footprint_rectangle` instead.\n",
      "  strel = {2: {'min': disk(1), 'max': square(3)}, 3: {'min': ball(1), 'max': cube(3)}}\n",
      "/Users/amin/Code/OpenPNM/.venv/lib/python3.13/site-packages/porespy/metrics/_funcs.py:65: FutureWarning: `square` is deprecated since version 0.25 and will be removed in version 0.27. Use `skimage.morphology.footprint_rectangle` instead.\n",
      "  strel = {2: {\"min\": disk(1), \"max\": square(3)}, 3: {\"min\": ball(1), \"max\": cube(3)}}\n",
      "/Users/amin/Code/OpenPNM/.venv/lib/python3.13/site-packages/porespy/metrics/_funcs.py:65: FutureWarning: `cube` is deprecated since version 0.25 and will be removed in version 0.27. Use `skimage.morphology.footprint_rectangle` instead.\n",
      "  strel = {2: {\"min\": disk(1), \"max\": square(3)}, 3: {\"min\": ball(1), \"max\": cube(3)}}\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"><span style=\"color: #7fbfbf; text-decoration-color: #7fbfbf\">[00:33:41] </span><span style=\"color: #800000; text-decoration-color: #800000; font-weight: bold\">ERROR   </span> Loaded version from importlib.metadata                                           <a href=\"file:///Users/amin/Code/OpenPNM/.venv/lib/python3.13/site-packages/porespy/__init__.py\" target=\"_blank\"><span style=\"color: #7f7f7f; text-decoration-color: #7f7f7f\">__init__.py</span></a><span style=\"color: #7f7f7f; text-decoration-color: #7f7f7f\">:</span><a href=\"file:///Users/amin/Code/OpenPNM/.venv/lib/python3.13/site-packages/porespy/__init__.py#57\" target=\"_blank\"><span style=\"color: #7f7f7f; text-decoration-color: #7f7f7f\">57</span></a>\n",
       "</pre>\n"
      ],
      "text/plain": [
       "\u001b[2;36m[00:33:41]\u001b[0m\u001b[2;36m \u001b[0m\u001b[1;31mERROR   \u001b[0m Loaded version from importlib.metadata                                           \u001b]8;id=172077;file:///Users/amin/Code/OpenPNM/.venv/lib/python3.13/site-packages/porespy/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=173415;file:///Users/amin/Code/OpenPNM/.venv/lib/python3.13/site-packages/porespy/__init__.py#57\u001b\\\u001b[2m57\u001b[0m\u001b]8;;\u001b\\\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import openpnm as op\n",
    "import porespy as ps\n",
    "import numpy as np\n",
    "from pathlib import Path"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5dceea09",
   "metadata": {},
   "source": [
    "Load a small image of Berea sandstone:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e83f28be",
   "metadata": {},
   "outputs": [],
   "source": [
    "pkg_root = Path(op.__file__).resolve().parents[2]\n",
    "path = pkg_root / \"tests/fixtures/berea_100_to_300.npz\"\n",
    "data = np.load(path.resolve())\n",
    "im = data['im']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "429e869f",
   "metadata": {},
   "source": [
    "Note meta data for this image, which we'll use to compare the network predictions too later:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0000735d",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = {\n",
    "    'shape': {\n",
    "        'x': im.shape[0],\n",
    "        'y': im.shape[1],\n",
    "        'z': im.shape[2],\n",
    "    },\n",
    "    'resolution': 5.345e-6,\n",
    "    'porosity': 19.6,\n",
    "    'permeability': {\n",
    "        'Kx': 1360,\n",
    "        'Ky': 1304,\n",
    "        'Kz': 1193,\n",
    "        'Kave': 1286,\n",
    "    },\n",
    "    'formation factor': {\n",
    "        'Fx': 23.12,\n",
    "        'Fy': 23.99,\n",
    "        'Fz': 25.22,\n",
    "        'Fave': 24.08,\n",
    "    },\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb2f0380",
   "metadata": {},
   "source": [
    "Perform extraction using `snow2` from `PoreSpy`.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b2a3d3d4",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/amin/Code/OpenPNM/.venv/lib/python3.13/site-packages/porespy/filters/_snows.py:446: FutureWarning: `cube` is deprecated since version 0.25 and will be removed in version 0.27. Use `skimage.morphology.footprint_rectangle` instead.\n",
      "  peaks_dil = spim.binary_dilation(input=peaks_i, structure=cube(3))\n",
      "/Users/amin/Code/OpenPNM/.venv/lib/python3.13/site-packages/porespy/filters/_snows.py:594: FutureWarning: `cube` is deprecated since version 0.25 and will be removed in version 0.27. Use `skimage.morphology.footprint_rectangle` instead.\n",
      "  labels, N = spim.label(peaks > 0, structure=cube(3))\n",
      "/Users/amin/Code/OpenPNM/.venv/lib/python3.13/site-packages/porespy/tools/_morphology.py:174: FutureWarning: `cube` is deprecated since version 0.25 and will be removed in version 0.27. Use `skimage.morphology.footprint_rectangle` instead.\n",
      "  strel = cube(w)\n"
     ]
    }
   ],
   "source": [
    "snow = ps.networks.snow2(\n",
    "    phases=im,\n",
    "    voxel_size=data['resolution'],\n",
    "    boundary_width=[3, 0, 0],\n",
    "    accuracy='standard',\n",
    "    parallel_kw=None,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fd628f0",
   "metadata": {},
   "source": [
    "For more details on the use of the `snow2` function see [this example](https://porespy.org/examples/networks/reference/snow2.html).  The `snow2` function returns a `Results` object with results attaches as attributes.  We can see these attributes if we print it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d6fbf707",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "Results of snow2 generated at Fri Dec  5 00:33:54 2025\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "network                   Dictionary with 26 items\n",
      "regions                   Array of size (206, 200, 200)\n",
      "phases                    Array of size (206, 200, 200)\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n"
     ]
    }
   ],
   "source": [
    "print(snow)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22e3da7e",
   "metadata": {},
   "source": [
    "We can open the network in OpenPNM easily as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "0263c5f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "pn = op.io.network_from_porespy(snow.network)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "473028d2",
   "metadata": {},
   "source": [
    "It's usually helpful to inspect the network by printing it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f1bd687f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "══════════════════════════════════════════════════════════════════════════════\n",
      "net : <openpnm.network.Network at 0x11af6f6b0>\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  #  Properties                                                   Valid Values\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  2  throat.conns                                                  2638 / 2638\n",
      "  3  pore.coords                                                   1422 / 1422\n",
      "  4  pore.region_label                                             1422 / 1422\n",
      "  5  pore.phase                                                    1422 / 1422\n",
      "  6  throat.phases                                                 2638 / 2638\n",
      "  7  pore.region_volume                                            1422 / 1422\n",
      "  8  pore.equivalent_diameter                                      1422 / 1422\n",
      "  9  pore.local_peak                                               1422 / 1422\n",
      " 10  pore.global_peak                                              1422 / 1422\n",
      " 11  pore.geometric_centroid                                       1422 / 1422\n",
      " 12  throat.global_peak                                            2638 / 2638\n",
      " 13  pore.inscribed_diameter                                       1422 / 1422\n",
      " 14  pore.extended_diameter                                        1422 / 1422\n",
      " 15  throat.inscribed_diameter                                     2638 / 2638\n",
      " 16  throat.total_length                                           2638 / 2638\n",
      " 17  throat.direct_length                                          2638 / 2638\n",
      " 18  throat.perimeter                                              2638 / 2638\n",
      " 19  pore.volume                                                   1422 / 1422\n",
      " 20  pore.surface_area                                             1422 / 1422\n",
      " 21  throat.cross_sectional_area                                   2638 / 2638\n",
      " 22  throat.equivalent_diameter                                    2638 / 2638\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  #  Labels                                                 Assigned Locations\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  2  pore.all                                                             1422\n",
      "  3  throat.all                                                           2638\n",
      "  4  pore.boundary                                                         182\n",
      "  5  pore.xmin                                                             104\n",
      "  6  pore.xmax                                                              78\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n"
     ]
    }
   ],
   "source": [
    "print(pn)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dcd1152d",
   "metadata": {},
   "source": [
    "When we update the SNOW algorithm to `snow2` we removed many of the opinionated decision that the original version made.  For instance, `snow2` does not pick which diameter should be the *definitive* one, so this decision needs to be made by the user once they can the network into OpenPNM. This is illustrated below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f1aa6c19",
   "metadata": {},
   "outputs": [],
   "source": [
    "pn['pore.diameter'] = pn['pore.equivalent_diameter']\n",
    "pn['throat.diameter'] = pn['throat.inscribed_diameter']\n",
    "pn['throat.spacing'] = pn['throat.total_length']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6cc07911",
   "metadata": {},
   "source": [
    "The user also needs to decide which 'shape' to assume for the pores and throats, which impacts how the transport conductance values are computed.  Here we use the `pyramids_and_cuboid` models:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "f32e2f79",
   "metadata": {},
   "outputs": [],
   "source": [
    "pn.add_model(propname='throat.hydraulic_size_factors',\n",
    "             model=op.models.geometry.hydraulic_size_factors.pyramids_and_cuboids)\n",
    "pn.add_model(propname='throat.diffusive_size_factors',\n",
    "             model=op.models.geometry.diffusive_size_factors.pyramids_and_cuboids)\n",
    "pn.regenerate_models()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20104b10",
   "metadata": {},
   "source": [
    "More information about the \"size factors\" can be found in [this example](https://openpnm.org/examples/reference/simulations/size_factors_and_transport_conductance.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "590326c1",
   "metadata": {},
   "source": [
    "One issue that can happen in extracted networks is that some pores can be isolated from the remainder of the network.  The problem with 'disconnected networks' is discussed in detail [here](https://openpnm.org/examples/reference/networks/managing_clustered_networks.html).  The end result is that we need to first check the network health, then deal with any problems:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "7868024f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "Key                                 Value\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "headless_throats                    []\n",
      "looped_throats                      []\n",
      "isolated_pores                      [456, 574, 1115, 1158, 1179]\n",
      "disconnected_pores                  [456, 574, 1115, 1158, 1179, 1204, 1346, 1347, 1348]\n",
      "duplicate_throats                   []\n",
      "bidirectional_throats               []\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n"
     ]
    }
   ],
   "source": [
    "h = op.utils.check_network_health(pn)\n",
    "print(h)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6247f065",
   "metadata": {},
   "source": [
    "The above output shows there are both 'single' isolated pores as well as 'groups' of pores that are disconnected.  These all need to be \"trimmed\".  We can use the list provided by `h['disconnected_pores']` directly in the `trim` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "e6975f8c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "Key                                 Value\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "headless_throats                    []\n",
      "looped_throats                      []\n",
      "isolated_pores                      []\n",
      "disconnected_pores                  []\n",
      "duplicate_throats                   []\n",
      "bidirectional_throats               []\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n"
     ]
    }
   ],
   "source": [
    "op.topotools.trim(network=pn, pores=h['disconnected_pores'])\n",
    "h = op.utils.check_network_health(pn)\n",
    "print(h)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bdb57f8c",
   "metadata": {},
   "source": [
    "Now we are ready to perform some simulations, so let's create a phase object to compute the thermophysical properties and transport conductances:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "ebe5084d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"><span style=\"color: #7fbfbf; text-decoration-color: #7fbfbf\">[00:33:54] </span><span style=\"color: #808000; text-decoration-color: #808000\">WARNING </span> throat.entry_pressure was not run since the following property is missing:       <a href=\"file:///Users/amin/Code/OpenPNM/src/openpnm/core/_models.py\" target=\"_blank\"><span style=\"color: #7f7f7f; text-decoration-color: #7f7f7f\">_models.py</span></a><span style=\"color: #7f7f7f; text-decoration-color: #7f7f7f\">:</span><a href=\"file:///Users/amin/Code/OpenPNM/src/openpnm/core/_models.py#480\" target=\"_blank\"><span style=\"color: #7f7f7f; text-decoration-color: #7f7f7f\">480</span></a>\n",
       "<span style=\"color: #7fbfbf; text-decoration-color: #7fbfbf\">           </span>         <span style=\"color: #008000; text-decoration-color: #008000\">'throat.surface_tension'</span>                                                         <span style=\"color: #7f7f7f; text-decoration-color: #7f7f7f\">              </span>\n",
       "</pre>\n"
      ],
      "text/plain": [
       "\u001b[2;36m[00:33:54]\u001b[0m\u001b[2;36m \u001b[0m\u001b[33mWARNING \u001b[0m throat.entry_pressure was not run since the following property is missing:       \u001b]8;id=970589;file:///Users/amin/Code/OpenPNM/src/openpnm/core/_models.py\u001b\\\u001b[2m_models.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=629820;file:///Users/amin/Code/OpenPNM/src/openpnm/core/_models.py#480\u001b\\\u001b[2m480\u001b[0m\u001b]8;;\u001b\\\n",
       "\u001b[2;36m           \u001b[0m         \u001b[32m'throat.surface_tension'\u001b[0m                                                         \u001b[2m              \u001b[0m\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gas = op.phase.Phase(network=pn)\n",
    "gas['pore.diffusivity'] = 1.0\n",
    "gas['pore.viscosity'] = 1.0\n",
    "gas.add_model_collection(op.models.collections.physics.basic)\n",
    "gas.regenerate_models()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "609a0594",
   "metadata": {},
   "source": [
    "Finally we can do a Fickian diffusion simulation to find the formation factor:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "1943e69d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Formation factor of the extracted network is 22.33394375429123\n",
      "The compares to a value of 23.12 from DNS\n"
     ]
    }
   ],
   "source": [
    "fd = op.algorithms.FickianDiffusion(network=pn, phase=gas)\n",
    "fd.set_value_BC(pores=pn.pores('xmin'), values=1.0)\n",
    "fd.set_value_BC(pores=pn.pores('xmax'), values=0.0)\n",
    "fd.run()\n",
    "dC = 1.0\n",
    "L = (data['shape']['x'] + 6)*data['resolution']\n",
    "A = data['shape']['y']*data['shape']['z']*data['resolution']**2\n",
    "Deff = fd.rate(pores=pn.pores('xmin'))*(L/A)/dC\n",
    "F = 1/Deff\n",
    "print(f\"The Formation factor of the extracted network is {F[0]}\")\n",
    "print(f\"The compares to a value of {data['formation factor']['Fx']} from DNS\")\n",
    "np.testing.assert_allclose(F, data['formation factor']['Fx'], rtol=0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89681eea",
   "metadata": {},
   "source": [
    "And a Stokes flow simulation to find Permeability coefficient:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "4a500b65",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Permeability coefficient is 1.4732149604234128 Darcy\n",
      "The compares to a value of 1.36 from DNS\n"
     ]
    }
   ],
   "source": [
    "sf = op.algorithms.StokesFlow(network=pn, phase=gas)\n",
    "sf.set_value_BC(pores=pn.pores('xmin'), values=1.0)\n",
    "sf.set_value_BC(pores=pn.pores('xmax'), values=0.0)\n",
    "sf.run()\n",
    "dP = 1.0\n",
    "L = (data['shape']['x'] + 6)*data['resolution']\n",
    "A = data['shape']['y']*data['shape']['z']*data['resolution']**2\n",
    "K = sf.rate(pores=pn.pores('xmin'))*(L/A)/dP*1e12\n",
    "print(f'Permeability coefficient is {K[0]} Darcy')\n",
    "print(f\"The compares to a value of {data['permeability']['Kx']/1000} from DNS\")\n",
    "np.testing.assert_allclose(K, data['permeability']['Kx']/1000, rtol=0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8597118",
   "metadata": {},
   "source": [
    "Both of the above simulations agree quite well with the known values for this sample.  This is not always the case because network extraction is not always perfect. One problem that can occur is that the pore sizes are much too small due to anisotropic materials.  In other cases the pores are too large and overlap each other too much.  Basically, the user really need to double, then triple check that their extraction is 'sane'.  It is almost mandatory to compare the extraction to some known values as we have done above.  It's also a good idea visualize the network in Paraview, as explained [here](https://openpnm.org/examples/tutorials/10_visualization_options.html), which can reveal problems.  "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "openpnm",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
